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1. INTRODUCTION 

Swarm intelligence algorithms are applied in almost every area of science and engineering, since in 
practice, they have shown their high efficiency in solving, first of all, complex real-life engineering 
optimization problems [1-3]. Complex optimization problems are understood as high-dimensionality 
problems with a large-scale complex topology of the solution search space, nonlinear criteria and constraints. 
Also, these problems may have more than one objective function, stochastic and dynamic properties. A 
scientist with an understanding of various swarm algorithms’ mechanisms, their weak and strong sides can 
create hybrid algorithms, combining, as a rule, two swarm algorithms or swarm and evolutionary algorithms. 
It allows achieving a synergistic effect, increasing the overall efficiency of the resulting hybrid algorithm [4]. 
There are several directions for hybridization. 

Each algorithm in the ensemble works according to its own rules. However, they use they use the 
general population of agents or exchange the best-found solutions, for example, a hybrid of the PSO and the 
genetic algorithm (GA) [5], PSO + differential evolution [6], PSO + grey wolf optimizer [7], PSO + 
biogeography-based optimization [8]. This approach does not require a change in the underlying algorithms. 
They work in parallel, exchanging data, or using a shared (public) population. As a result, the complexity of 
the algorithms does not increase; the basic algorithms can be easily replaced. Simultaneously, in such an 
approach, the synergistic effect may be small because of the weak integration of algorithms into the system. 
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Adding a rule or stage of updating the positions of agents from another algorithm to the main 
algorithm. In [9], a mutation stage was added to the classical PSO algorithm at the end of each iteration, 
similar to that in GA; similarly, the mutation was added to the grey wolf optimizer in [10]; the study [11] also 
added a mutation stage to the PSO algorithm, but according to the principles of the differential evolution. 
Papers [12-13] describe PSO with an additional allowance for particle acceleration, calculated using formulas 
from the gravitational search algorithm. In this approach, a new algorithm is created, often using the 
principles of both swarm intelligence and evolutionary computation. 

Hybrid algorithms with a higher degree of integration of the rules and principles of operation of two 
algorithms are also possible: FireFly + PSO algorithm [14], PSO + artificial bee colony algorithm [15], PSO 
+ ant colony optimization (ACO) [16], PSO + wolf pack algorithm [17]. This method allows us to create an 
algorithm that considers problem’s features being solved and compensates for the shortcomings of one 
algorithm by adding properties of another to it. Nevertheless, this significantly complicates the created 
algorithm for solving the problem and selecting the heuristic parameters’ algorithm. In other words, the 
complexity of the hybrid algorithm exceeds the sum of the complexities of the parts that leave it. 

Using meta-optimization when one algorithm performs the selection of parameters of another 
algorithm. For example, in [18], local unimodal sampling was used to select the values of the PSO algorithm 
parameters; in the article [19], the PSO parameters are selected using GA; in [20], the ACO parameters are 
selected using the bacterial foraging algorithm. This approach requires a very large expenditure of 
computational resources. 

At each iteration of the algorithm, all agents' positions in the population or only the best agents are 
used as an initial approximation for the local search algorithm. In [21], the unstringing and stringing operator 
was used for the ACO with the local search; the use of various local search algorithms in PSO is discussed in 
[22]; FireFly + gradient descent is used in [23]; The PSO + lin-kernighan algorithm was considered in [24]. 

Swarm and evolutionary metaheuristic algorithms themselves are quite complex in understanding 
their work and fine-tuning to the features of the tasks being solved. Hybrids from different metaheuristic 
algorithms turn out to be even more complex and less universal. However, as shown in [17], the 
metaheuristic algorithm's efficiency can sometimes be increased on the contrary, by simplifying it. Perhaps 
that is why there are many works in which complex hybrid algorithms are successfully applied to solve a 
specific problem. Still, none of them has become as widespread as less complicated ACO and PSO among SI 
algorithms [4] or GA among evolutionary ones. 

Our approach to solving optimization problems aims to increase SI algorithms’ efficiency and 
versatility without increasing their complexity. We apply a combination of a universal, problem-free SI 
algorithm with simple deterministic domain-specific heuristic algorithms or, in other words, greedy 
heuristics. A hybrid of a stochastic universal SI algorithm and a deterministic simple domain-specific 
algorithm achieves a greater synergy effect than combining two swarm or evolutionary algorithms. 


2. GENERALIZED SCHEME OF SI ALGORITHM AND APPLICATION APPROACH 

To implement the proposed approach, it is necessary to make sure that it does not depend on the 
choice of a SI algorithm and a domain-specific heuristic algorithm. In this case, it will be quite versatile and 
flexible and will not require any modifications to the algorithms; only their interaction features will change. 
Therefore, it is necessary to define a generalized SI algorithm model so that various SI algorithms can be 
easily applied since it is impossible to say in advance which one will be better in a given task. Next, we need 
to define a general universal scheme of interaction between a SI algorithm and a domain-specific heuristic 
without strong dependencies. 


2.1. Generalized SI algorithm model 

An SI algorithm is specified by data structures and operations on them, like many other algorithms. 
SI uses the swarm, i.e., a population of agents as a basic data structure and rules of moving agents as an 
algorithm for transforming the data. A distinctive feature of the SI algorithms is using an indirect information 
exchange between agents. For example, the PSO uses the best-found position in the search space [25]; the ant 
colony optimization uses a pheromone graph [26]; the monkey search algorithm uses a weighted center of 
agents’ positions [27]. It is possible to formulate basic parts of an SI algorithm: 
— A set of agents S. 
— An object for the indirect interaction M. 
— An algorithm of agents’ moving and interaction with optimization problem A. 
— Heuristic parameters P. 
—  Aninput to obtain data from the optimization problem being solved J. 
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— An output to send solutions of the optimization problem and the final algorithm result O. 

Any SI algorithm works according to the following scheme. 

— Generation of the initial population. The swarm agents distribute randomly in the solution searching 
space. 

— Calculation of fitness-functions. Each agent receives the values of the optimization problem criterion 
(criteria). If the agent has a better fitness than all previous finesses of all agents, this solution is saved as 
the best current solution. 

— Agents’ movement. The agents move in the search space using indirect interactions. 

— If the stop condition is satisfied, the algorithm needs to be finished or otherwise needs to be passed to 
Step 2. The saved best solution is the algorithm result. 


2.2. Interaction between a SI algorithm and a domain-specific algorithm 

The following interface and interaction steps were proposed for the interaction of an SI algorithm 
and a domain-specific heuristic algorithm. 

— The SI algorithm receives a dimension of the search space of an optimization task via input J. It is used 
as the length of the vector X. 

— The SI algorithm generates the variants of the problem solution by agents’ movements. A count of 
variants is a count of agents. 

— Each position is sent to a domain-specific heuristic algorithm of an optimized system via output Opos. 

— The domain-specific heuristic algorithm encodes the position and uses it for solving an optimization 
problem and calculating the criterion values. Thus, the SI algorithm works as a meta-optimizer for the 
domain-specific algorithm. 

— All obtained criterion values are sent to the SI via /. 

— The SI algorithm receives criterion values and performs agents’ movement, i.e., we go to the second 
step of this algorithm. 

When a stop-condition is satisfied, then the best agent position is transmitted to output O. The 
simplest stopping condition is the execution of a given number of iterations. 

The interaction’s listed steps provide the independence of used algorithms. The approach proposed 
allows us to implement different SI algorithms quickly. As far as changes in the optimization problem do not 
cause the necessity of changing anything in the SI algorithm implementation, and vice versa. The absence of 
close relationships makes it possible to combine any SI algorithm with domain-specific heuristics. 


2.3. Mapping of agents’ positions 

In the proposed approach, we use the search space constraints from 0 to 1 for each dimension. The 
introduced restrictions are necessary for successful adaptation and simple integration of the algorithm for 
solving various optimization problems. The fact is that some SI parameters have a different effect depending 
on the scale of optimization variables. For example, in the firefly optimization algorithm [28], the distance 
between agents affects agents’ attraction nonlinearly. Therefore, if a range of search space is from 0 to 100, 
the distance from the middle (50) to the high edge (100) will be 75. While if this distance is scaled from 0 to 
1, then the same distance will be equal to 0.75, and the degree of attraction between the agents would be 
completely different. At the same time, this distance is equivalent to the optimization task. Also, the radius of 
the neighborhood region width of the bee’s algorithm [29]. The larger the distance between the boundaries of 
a search space direction, the greater the neighborhood region width should be. It would be necessary for a 
non-homogeneous search space to introduce different values of the parameters for different directions. 

As a result, the algorithm parameters would have to be changed even if a simple change of 
measurement units in the task (hours to seconds, and kilometers to feet). We suggest specifying the search 
space bounded from 0.0 to 1.0 in all directions into the algorithm to avoid this. And mapping of the agents’ 
coordinates from the algorithm into the values of the optimized variables. In the simplest variant, this can be 
done as follows. We denote the coordinate vector of an agent X = {x1, x2, ..., x.}, and the vector of optimized 
variables Q = {q1, q2, .... gz}. Under the restrictions a; < gi < bi, we obtain a map of the form 
qi = xi(bi — ai) +a,i=1,...,L. 

A more sophisticated variant is shown in sections 3-4. The proposed parameter mapping is 
necessary to eliminate unnecessary relationships between the optimization problem model and the 
optimization algorithms. 
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3. RESULTS AND DISCUSSION 
3.1. Application example. job-shop scheduling problem 

This subsection introduces its own symbolic notation, which should not be confused with the 
notation in the descriptions of the SI algorithms (Section II). As mentioned above, a model of an optimization 
problem and the SI algorithms are completely independent. The job-shop scheduling problem is among the 
most challenging combinatorial optimization problems. Scheduling tasks may be characterized as one of the 
most significant optimization problems since plans and schedules need to be arranged in all fields. The 
problem consists of scheduling a set of jobs N on a set of machines M to minimize the processing time named 
as the makespan. It is the time needed to process all jobs. Each job includes a number of operations and has a 
predefined operation order. The order defines a specific machine and time interval for each operation and 
sequence of operation. Each machine can process only one operation at a time. It needs to find the shortest 
(quickest) schedule as a distribution of the operations to time intervals on the machines. This problem 
belongs to the NP-hard class [30]. Let: N = {1, ..., 1} is the set of jobs; M= {1, ..., m} is the set of machines; 
O = {00, ..., o1+1} denote the set of the all operations, 0 and L+1 are fictions operations: start and finish, L is 
the total number of operations of all jobs [31-32]. 

Each operation o € O has its processing time p(o). And A denotes the set of ordered pairs of 
operations constrained by the precedence relations; FE; denotes the set of all pairs of operations which require 
machine M,. Next, let t(0) is the start time of the operation 0, t(09) = 0. The problem is to find a schedule: 


max(t(0) + p(o)), 0 € O; (1) 
t(0) > 0, Vo € O; (2) 
t(v) — (0) = po), W(o, v) € A; (3) 
t(w) — t(o) 2 p(o) V t(o) — t(w) 2 pw), (4) 


Viw,o) € Ek k=1,..., m. 
The value of max(t(o) + p(o)) is called the makespan. 


3.1.1. SI Algorithm with the greedy heuristic for job-shop scheduling problem 

The following algorithm can be used as a simple greedy heuristic for solving the problem. 
Input: O, A, L, Ex, m,n 
1. Ve {} 

2. Qi: 1-1 <— sort(O) in descending order and preserving the order given in Ex, k = 1,...,m 
3. 20 <— 00, QL<— OL 
4. (Qo) — 0, ¥ — {a0} 
5.Fori=1...L: 
5.1. place the operation Q;: 
t(Q;) — min(t(Q;)) under conditions: 
t(v) — (Qi) = p(Qi), Vv) € ¥, W(Qi, v) € A; 
t(w) — t(Q)) = p(Qi) V (Qi) — tw) = pw), Vw) € Y% Vw, Q) € Ei, k= 1, ...,m 
5.2. ~— PuU{Q;} 
Output: Y (vector of distributed operations) 

As a result, for each operation o € O, the start time f(o) will be determined, and the value of the 
optimality criterion will be equal to 7(0z). 

This greedy heuristic allows us to schedule activities to prioritizing more extended activities. The 
longer the execution time of the operation p(o), the more difficult it is to find a free machine for it; therefore, 
priority is given to more extended operations. The advantage of the algorithm is its simplicity and the ability 
to add additional restrictions on the relationship of operations to it. However, such an algorithm can give 
solutions that are far from optimal. 

The algorithm's efficiency in terms of the criterion of the problem (minimization of the makespan) 
can be significantly increased due to the intelligent prioritization of operations. In this case, the same 
heuristic scheme can be applied only to order the operations not by their duration but with priorities, which 
can be determined using SI algorithms. For applying an SI algorithm to the job-shop scheduling problem, an 
agent position (vector X) is needed to map into a possible schedule. The vector X must contain as many 
elements as total operations in all jobs (L). The mapping process consists of several steps. 


Int J Artif Intell, Vol. 10, No. 1, March 2021: 157 — 165 


Int J Artif Intell ISSN: 2252-8938 i) 161 


Input: O, X 

1. OX — concatenate_element_by_element(O, X) 
2. OX — sort(OX) ascending x. 

3. O* — OX): 1,1 (remove x from pairs {o, x} 


4.Fori=1...L: 
replace the operation 0; by the operation 0, that should be on the correspondence order for this job. 
Output: O* 


Consider an example of 3 works, each of which consists of 3 operations. 

J = {01, 02, 03}; J2 = {04, 05, 06}; J3 = (07, 08, Oo}. 

Let X = {0.11, 0.72, 0.57, 0.92, 0.43, 0.21, 0.12, 0.40, 0.3}. 

Then O= {i, Jo, J3} = {01, 02, 03, 04, 05, 06, O7, O8, 09}. After sorting: 

OX = {{04, 0.92}, {02, 0.72}, {03, 0.57}, {05, 0.43}, {08, 0.4}, {09, 0.3}, {06, 0.21}, (07, 0.12}, {o1, O.11}}. 
The operation order obtained cannot be used as a schedule since it is possible to violate the sequence of 
operations within one job. 

Finally, taking into account the restoration of order within the works, we obtain O* = {00, 04, 01, 02, 
05, 07, 08, 09, 06, 03, O10}, Oo and ojo are fictions start and finish operations. 

The vector O* does not exactly mean the order of the starting of the stages. It means that operations 
are extracted from O in the order in which they are located in O*. For the extracted operation, the moment of 
the fastest possible start is determined. Thus, the vector O° does not specify the order of processing 
operations, but the order of determining the start time for them. Since the time costs for performing each 
operation are assumed to be constant, the vector O* uniquely determines the final schedule. Therefore, the 
vector X is uniquely translated into the required discrete schedule. The advantage of the above scheme is 
obtaining an allowable schedule for any values of the vector X. 


3.1.2. Experimental results 

The PSO algorithm was used to solve the set of job-shop problem instances [33] denoted as LAOI- 
LA21 of various sizes provided by OR-Library [34] (LA because the author is S. Lawrence). The sizes of 
instances are from 10x5 to 15x10 (nxm). Table 1 shows the results (average and the best solutions from ten 
runs of the algorithm). Also, Table | list the best well-known results from the literature [34]. The values of 
the PSO parameters: population size 50, a = 1.76, a2 = 1.38, w = 0.73, 6 = 0.28 [18]. 

The summary deviation between the best and well-known results is 135 hours (0.71%) using the 
PSO. The results are close to the best (deviation < 1%) for most tasks (81%). Thus, the PSO algorithm allows 
solving job-shop scheduling problems very close to the best-known results. Moreover, the PSO algorithm is 
high-speed, since only 100 PSO iterations were used. Increasing the iteration number can give better results 
for a longer time. The largest time of solving one instance was about 0.5 seconds using 2.4 GHz Intel CPU 
i7. The minimum time was 0.08 seconds. Such variation is explained by the different dimensions of the 
instances. 


Table 1. Results of the PSO with the greedy heuristic in job-shop scheduling problem instances 


Instance Avg Best Best known 
LAOI 690.8 666 666 
LAO2 693.35 658 655 
LAO3 541.49 597 597 
LA04 613.79 590 590 
LAOS 593 593 593 
LA06 926.03 926 926 
LAO7 902.8 890 890 
LA08 882.09 863 863 
LAO9 955.34 951 951 
LA1O 958 958 958 
LAI1 1222.05 1222 1222 
LA12 1043.59 1039 1039 
LA13 1161.59 1150 1150 
LA14 1292 1292 1292 
LAI5 1283.31 1207 1207 
LA16 1011.11 956 945 
LA17 819.67 784 784 
LA18 930.19 859 848 
LAI9 865 865 842 
LA20 907 907 902 
LA21 1128 1128 1046 
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3.2. Application example. exploring potential energy surface of nanocluster isomers 

This subsection discusses two related problems. The first problem is the local minima search on the 
potential energy surface, which are isomers for a cluster of a given composition. The second problem is 
calculating the «reaction path» as the energetically most favorable path of the phase point movement between 
two known local minima. The lower the minimum, the more stable corresponding configuration of atomic 
nuclei. There are specialized software solutions in the field of bioinformatics for exploring ligand 
binding/unbinding pathway [35]. Isomers are two (or more) molecules (nanoclusters) with the same 
molecular formula, i.e., atomic composition. Isomers are divided into two main categories. Structural isomers 
(or constitutional isomer) have the same formula, but the atoms are bonded together in different orders. 
Stereoisomers have the same connectivity but the different spatial arrangements [36]. Isomers can have very 
different physical properties, such as boiling point, melting point, and chemical reactivity. In metal systems, 
the properties of stereoisomers, as a rule, do not differ. 

The act of a structural phase transition in a nanocluster is a mutual geometric rearrangement of 
atoms occurring in time [37]. Accordingly, the cluster's energy and its physicochemical properties 
substantially depend on the geometric arrangement of the atoms. Even insignificant changes in the spatial 
structure of a cluster can lead to rather large changes in its total energy. The initial geometry's suitable choice 
is crucial for obtaining reliable, calculated values of a cluster's physicochemical characteristics in a stationary 
(equilibrium) state. For a successfully calculated geometry, changes in energy values during the optimization 
process should be insignificant (less than a given value). Also, for an optimized geometry, the first energy 
derivatives (gradient) with respect to all geometric distortions should be close to zero. It suggests that we are 
on a flat energy surface. When considering the broader features of the potential energy surface (PES), 
including topology and topography, it has become usual to refer to the PES as the «potential energy 
landscape» [38]. 

Energy surfaces and landscapes hold the key to understanding a wide range of molecular 
phenomena. To construct the PES, it is convenient to use the so-called internal coordinates, which include 
bond lengths or other interatomic distances, bond and dihedral angles of a cluster/molecule. The number of 
3n-6 (where n is the number of atomic nuclei) of independent coordinates should include those structural 
parameters that change most dramatically during the transformation under investigation. The paper [39] 
describes a potential energy surface in terms of local minima and the transition states that connect them 
provides a conceptual and computational framework for understanding and predicting observable properties. 

A potential barrier to structural transformation is the energy of the transition state (saddle point) 
relative to the initial minimum of the potential energy surface. It is the potential barrier that determines the 
rate of the process of structural transformation. In the Arrhenius equation, these are the activation energies 
Ex: k = koexp(-Ea/ kpT); Ea is the energy activation, ko is the pre-exponential factor for the reaction, kg is the 
Boltzmann constant, T is the absolute temperature (in Kelvins). 

To construct a minimum energy path between two known constitutional isomers, an efficient 
solution to the assignment problem in a bipartite graph is required. Each atom of isomer A (reagent) is 
associated with one and only one atom of isomer B (product). The problem can be solved as a global 
optimization problem. Applying traditional numerical methods is impractical because they need huge 
amounts of computational resources like time and memory [40]. 


3.2.1. Greedy heuristic 

Thus, the search for isomers corresponding to local minima and the construction of the minimum 
energy path between them are the basic stages of the PES construction process. It is necessary to find the 
permutation function P, which would associate with each atom a; from A (i= 1, ..., 2) one and only one atom 
by) from B, so that root mead squared distance, RMSD(f)=(2/"z1|ai— b;'|?)'!? min, b;*=bp.y 

Additional restrictions are imposed on collisions of the transition paths of atoms. The minimum 
distance between the centers of the vectors connecting the atoms aj and bpy (i = 1, ..., 2) must be no less than 
Cmin = Cais MIn(min_disance_A, min_disance_B), where Cais is an empirical coefficient (~0.8), min_disance_X 
— the minimum distance among all pairs of atoms of isomer X (A or B). This limitation, high dimensionality 
and the presence of several types of atoms do not allow the existing methods for solving assignment 
problems without their significant modification. 

We combine a straightforward greedy heuristic algorithm and the universal PSO algorithms. The 
heuristic used involves sequentially sorting the atoms of isomer B and matching each of them with the 
nearest atom from A that has not yet been matched. This heuristic can be written: 

Input: A, B, n, L (L list of numbers from 1 to n without repetitions) 
1. tabu — {}. 
2.Fori=1...nj=1...0: 

i — |ai— bj? (for atoms of different sorts, the distance is infinite) 
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3. Fori=1...1: 

3.1. j<— argmin(D,;;), j € tabu 
3.2. P(Li) —j 

3.3. tabu <— tabu U {j}. 
Output: P 

Step 3.2 means that L; atom of B goes over to the position of the j-th atom of A. 

The resulting table function P is the desired permutation. The advantage of such heuristic is a very 
high speed of operation so the solution is obtained in one pass through the atoms of isomer B. If the phase 
structure of two isomers is close and isomer B is not rotated relative to isomer A, this solution method allows 
to find the optimal solution regardless of the number of atoms and number sorts of atoms. But it is not 
suitable for more complex problems since at the very first iterations of the work, it can match atoms of 
isomer B to such atoms of isomer A, which, if optimally solved, should be matched with other atoms of 
isomer B. 


3.2.2. SI algorithm with the greedy heuristic 

The resulting algorithm can be written: 
Input: A, B,n, m, pso_iters 
1. Xn — {x1, X2, ..., Xn}, Xmerandom(0, 1), m= 1, ..., |S], A=1, ...,7 
2. For i= 1... pso_iters: 

2.1. For m=1, ..., |S]. 

2.1.1. L — create_order_list {Xin 1, Xm2, .--) Xmn}- 
2.1.2. P — greedy_heuristic(A, B, n, L) 

2.1.3. fitness — RMSD(A, P(B)) 

2.1.4. If fitness < fitness_min 

fitness_min < fitness 

Poest <P 

2.2. Particles’ movement 

Output: Pres: 

The developed method has been tested on some isomer configurations from 7 to 39 atoms. Table 2 
shows the values of the RMSD obtained by the greedy heuristic only and by the PSO with the greedy 
heuristic. The mathematical modeling confirmed that the combination of greedy heuristics based on the 
physical properties of the problem and SI algorithms significantly improved the results regarding applying 
these methods separately. 


Table 2. Results of finding the shortest distance between isomers of the bimetallic cluster 
Instance Greedy heuristic, RMSD, Angstrom _ PSO with the greedy heuristic, RMSD, Angstrom 


Ag7 9,78 9,02 
Aul9 24,52 20,80 
Au20 no valid solution found 26,88 

Aul2Ag12 no valid solution found 98,41 
Cul5Ag15 63,85 61,72 
Col6Cul6 no valid solution found 67,61 
Au6Ag33 no valid solution found 48,17 


4. CONCLUSION 

In this article, a hybrid approach is proposed utilizing the strengths of SI and domain-specific 
heuristic algorithms. The main idea behind developing is: (a) to help the SI algorithm to take into account the 
specifics of the solved problem thanks to the domain-specific algorithm and (b) to improve the efficiency of 
the domain-specific algorithm by introducing stochastic and self-organized properties from the SI algorithm. 
Different SI algorithms can be easily applied due to SI algorithm’s low coupling and the domain-specific 
heuristic algorithm. The approach’s application is demonstrated on the job-shop scheduling problem and the 
problem of construction of the minimum energy path between two isomers. The experiments confirmed that 
the proposed approach allows the use of SI algorithms as a meta-optimizer that increases domain-specific 
heuristic algorithms’ efficiency. 
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